A numerical approach to the ground and excited states of a Bose-Einstein condensed 

gas confined in a completely anisotropic trap 
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The ground and excited states of a weakly interacting and dilute Bose-Einstein condensed gas, 
confined in a completely anisotropic harmonic oscillator potential, are determined at zero temper- 
ature within the Bogoliubov approximation. The numerical calculations employ a computationally 
efficient procedure based on a discrete variable representation (DVR) of the Hamiltonian. The 
DVR is efficient for problems where the interaction potential may be expressed as a local function 
of interparticle coordinates. In order to address condensates that are both very large (millions of 
atoms) and fully anisotropic, the ground state is found using a self-consistent field approach. Ex- 
perience has demonstrated, however, that standard iterative techniques applied to the solution of 
the non- linear partial differential equation for the condensate are non-convergent. This limitation 
is overcome using the method of direct inversion in the iterated subspace (DIIS). In addition, the 
sparse structure of the DVR enables the efficient application of iterative techniques such as the 
Davidson and/or Lanczos methods, to extract the eigenvalues of physical interest. The results are 
compared with recent experimental data obtained for Bose-Einstein condensed alkali metal vapors 
confined in magnetic traps. 
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I. INTRODUCTION 

The experimental achievement of Bose-Einstein con- 
densation (BEC) in dilute alkali-metal gases confined 
in magnetic traps Jj]-^| has generated tremendous in- 
terest in the behavior of the inhomogeneous, weakly 
interacting, and dilute Bose gas. At low tempera- 
tures, the confined Bose gases have been shown to 
be well-described by mean-field linear-response theories 
(MFLRT) based on the Bogoliubov |^,^] approximation, 
which assumes that the number of condensate atoms 
Nq is a substantial fraction of the total number [f^-fTs]. 
The finite-temperature extensions of this theory, such as 
the Hartree-Fock-Bogoliubov and Popov approximations 
Jj^ , fl5|], a lso have been successfully applied to these sys- 
tcms jL6|-p2[, though the microscopic basis for these ap- 
proaches and their agreement with experiment remain 
somewhat uncertain [p2[-p4|. 

The usual approach taken within MFLRT is to 
first solve the (time-independent) non-linear Schrodinger 
equation for a given number of atoms in the conden- 
sate; the resulting wavefunction and chemical potential 
are then used in the linear(ized) equations for the quasi- 
particle excitations in order to obtain the eigenmodes of 
the system. At finite temperatures, the significant de- 
pletion of the condensate requires that this procedure be 
iterated to self-consistency. The magnitude of the non- 
linear, or self-energy, term appearing in the Schrodinger 
equation for the ground state is proportional to the con- 
densate density. Thus, a fully three-dimensional solution 
of the non-separable equation for the condensate presents 



a number of numerical challenges, particularly for large 
values of No- 

The magnetic trap used to confine these gases is well 
approximated by a harmonic potential, in which the fre- 
quencies (uj x , Lu y , lu z ) are generally incommensurate. All 
of the published experimental realizations to date have 
employed traps with either spherical or cylindrical sym- 
metry, which are more conducive to numerical study. 
The number of experimental groups which have obtained 
BEC of confined alkali-metal gases is steadily rising, 
however p5[ ]. In order to consider a range of possible 
trap geometries, as well as to permit future investiga- 
tions of any time-dependent properties, it is necessary to 
generalize the numerical calculations to allow for com- 
plete anisotropy. The purpose of the present paper is to 
demonstrate that the numerical obstacles may be over- 
come and to outline a robust and computationally effi- 
cient procedure within the MFLRT. 

The approach taken by us is based on a discrete vari- 
able approximation (DVR) (2^] to the equations govern- 
ing the statics of the condensate. The DVR has several 
advantages over the methods used by others for the cur- 
rent problem. Most standard techniques discretize the 
solution to the partial differential equations (PDE) for 
the condensate and excitations either in physical (grid) 
space or in function (basis set) space. Grid-based meth- 
ods |27[] have the advantage that all local interactions 
between particles have a local representation while the 
kinetic energy takes on a sparse and typically banded 
structure. This sparsity is a consequence of (typically) 
low-order finite difference approximations to the deriva- 
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tives in the PDE, and is crucial for the implementation 
of iterative techniques 29 developed for large, sparse 
linear systems 30 1. The disadvantage of grid methods 



condensate and the Bogoliubov |9j6| equations for the 
excitations: 



is that it is often difficult to approximate the derivatives 
that appear in the PDE's to sufficient accuracy with- 
out resorting to high-order differences or very small step 
sizes. While an accurate representation of the kinetic en- 
ergy is arguably of minor importance for the ground state 
when No is large, it is crucial for the subsequent determi- 
nation of excitations. Expansions in function space (often 
called spectral and/or pseudospectral methods |31|) are 
typically superior since it is possible to "analytically" dif- 
ferentiate the functions without further approximation. 
These methods, however, require the evaluation of po- 
tential energy matrix elements by quadrature, leading to 
a non-diagonal and dense matrix representation of local 
operators. 

The DVR exploits the dual relationship between cer- 
tain orthogonal polynomials (such as the classical orthog- 
onal polynomials) and the points and weights of a Gauss 
quadrature. By using an appropriate set of polynomi- 
als as a basis for expanding the solution to the PDE, 
it is possible to maintain almost all of the advantages 
of a grid-based approach as well as the global conver- 
gence of a finite basis set. In the DVR, any local po- 
tential energy operator is diagonal and therefore easy to 
compute. The multidimensional kinetic energy also has 
a sparse representation because it is a commuting sum 
of one-dimensional operators. The sparseness is not as 
structured as is the case for more traditional approaches, 
but the expense of generating the kinetic energy matrix 
elements is mitigated by the iterative method used to 
solve the PDE: these terms need only be evaluated once. 
It is also worth noting that the expansion coefficients of 
the solutions to the PDE are trivially related to the val- 
ues of the solution on the quadrature points which serve 
as the physical grid. Perhaps the most important point 
to stress, though, is that the solutions of the PDE scale 
very efficiently with basis set size. This will become very 
important for three-dimensional problems where the size 
of the various matrices could easily be as large as 100 000 
by 100 000. 

The usual starting point for the theoretical treatment 
of the inhomogeneous, weakly interacting, dilute Bose 
gas at zero temperature is the time-independent MFLRT 
in the Bogoliubov approximation. Since the derivation of 
the resulting equations may be found in many places Q , 
only the central results are presented here. In the absence 
of inhomogeneities other than the confining potential, all 
wavefunctions may be assumed to be real. In the weakly 
interacting and low-density limit relevant to the exper- 
iments of interest, the interatomic interactions may be 
approximated by a two-body delta-function contact po- 
tential with a coupling constant g. Minimizing the grand 
canonical potential for interacting bosons and then lin- 
earizing the resultant equations for the amplitude of the 
condensate ip(r) and the excitations, u(r) and -u(r), leads 
respectively to the Gross-Pitaevskii (GP) p3|,|9| for the 



Lip{r) = /#(r); 



(1) 



(l - /j, + Vn) u n (r) - VaV n (r) = e n u n (r); (2) 
L — n + Vn) v n (r) — Vnu n (r) = -e n v n (r), (3) 



where the non-linear Schrodinger operator is written as 



2M 



V 2 + Vtrap(r) + Va. 



(4) 



The non-linearity of the GP equation is due to the mean- 
field, or Hartree, potential Vh = 5l?M r )| 2 ; a ^ l° n g wave- 
lengths, the coupling constant g « 4irh 2 a/M, where a 
and M are respectively the s-wave scattering length and 
the mass of the atom. Only the case a > will be consid- 
ered here. The chemical potential, fixes the number of 
atoms Nq in the condensate (the contribution from the 
excitations is ignored in the Bogoliubov approximation) , 
and e n are the collective excitations of the system. Note 
that the condensate is the zero-energy solution of the Bo- 
goliubov equations, ip(r) = u (r) = t>o(r); the excitation 
frequencies are measured with respect to the ground state 
energy fi. 

The confinement due to the magnetic trap is well- 
described by a completely anisotropic harmonic oscillator 
potential: 
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(5) 



The analysis of the GP and Bogoliubov equations is con- 
siderably simplified by re-expressing the trap frequen- 
cies u and coordinates in terms of the scaled variables 
(w x ,u) y ,u) z ) = tL>o(l,a, /?) and (x,y,z) -> d {x,y,z), 
where do = y/h/ Mujq is the characteristic oscillator 
length in the i-direction. As a result, the coupling con- 
stant becomes g — » 47r?7o, where we define r]o = Noa/d®. 
The Schrodinger operator is then 

L = -iv 2 + X - (x 2 + a 2 y 2 + (3 2 z 2 ) + A7r Vo ^ 2 (r). (6) 

All energies (including fj,,e n ) are now given in trap 
units fkOQ. With these choices, the condensate and 
excited-state wavefunctions are normalized to unity in 
the rescaled (dimensionless) coordinates: 



(7) 

(8) 



dr V> 2 (r) = 1; 
dr [u 2 n (r) - « 2 (r)] = 1. 



In the next section we discuss the numerical methods 
used to solve these equations. 



2 



II. NUMERICAL METHODS 

The discrete variable representation is particularly use- 
ful for the problem of trapped Bose-condensed atoms. As 
discussed in the Introduction, all potential energy oper- 
ators in the Hamiltonian are local and therefore have a 
diagonal matrix representation, while the kinetic energy 
still has a fairly simple and sparse structure. Since only 
the low-lying modes are relevant to the static proper- 
ties of the weakly interacting Bose gas at low tempera- 
tures, the sparse Hamiltonian matrices are ideally suited 
to iterative techniques for the determination of eigenval- 
ues. Since these methods are dominated by matrix- vector 
multiplies, a structured and sparse matrix offers consid- 
erable computational savings. 

An important feature of the GP equation is that the 
interaction potential depends non-linearly on the solu- 
tion. Thus, for a given Nq, the lowest eigenvalue fi and 
eigenvector ^>(r) must be found self-consistently. In prac- 
tical terms, a robust iterative method is required. Indeed, 
as discussed below, no solution to the GP equation may 
be found simply by straight iteration when the universal 
scaling parameter rj = aPr/o ^ 10. Since a/do is typi- 
cally of order 10 -3 — 10~ 2 , only a few thousand atoms 
in the trap (several orders of magnitude fewer than are 
experimentally relevant) would be sufficient to prevent 
the solution to the GP equation by iterative methods. 
We have developed a variant of a well-known technique, 
called direct inversion in the iterative subspace (DIIS) 
p5| to render convergent the self-consistency process for 
the large Nq of practical interest. DIIS complements al- 
ternative numerical approaches, such as the method of 
steepest descents (imaginary-time propagation), which 
have been successfully applied to the trapped Bose con- 
densates [B6^ 



A. DVR Techniques 

The basic ideas of a discrete variable representation are 
quite old. The earliest applications |39| were designed to 
simplify the calculation of certain classes of matrix ele- 
ments that appeared in finite basis set, yariational cal- 
culations. Subsequent authors [§(|[io] [l3| began to use 
the DVR more directly and, in some instances, to view 
it as the primary representation for the problem under 
consideration. The viewpoint of the current authors is 
that the DVR follows naturally from a particular choice 
of a finite basis set, one which is mathematically linked 
to Gauss quadratures. 

Let us consider a basis of functions, {<fi n (x), n — 1, N}, 
perhaps satisfying some set of boundary conditions over 
a finite or infinite interval, which is complete enough for 
expanding any unknown function over that interval to 
sufficient accuracy. For simplicity we assume that these 
functions form an orthonormal basis for the space. What 



we seek is a complementary set of "coordinate eigenfunc- 
tions", {v,i(x), i = 1, N}, and a generalized quadrature 
consisting of roots and weights {xk] Wk,k — 1, N}, such 
that 



Ui{x k ) = S ik . 



(9) 



We now expand the unknown Ui(x) in the set of functions 



N 



H{x) = ^2 <t>n{%)(<fin | Ui) 



(10) 



and assume it is possible to evaluate the overlap integral, 
((/) n | Ui), sufficiently accurately using the generalized 
quadrature rule so that 



b n | Ui) = Wi<fi n (Xi)\ 



N 



Ui(x) =W,^ (fi n (x)(j) n (xi). 



(11) 

(12) 



The result (|llj) follows directly from Eq. (|9|) and the in- 
tegration rule for a Gauss quadrature: 

,& 

(f\g)= dx w(x)f(x)g(x) = ^2 w kf(xk)g(x k ), 
Ja k 

(13) 

where w(x) is a non-negative weight function. Moreover, 
the delta-function property ^ of the coordinate eigen- 
functions is exactly satisfied by Eq. ( |l2| ) at the quadra- 
ture points Xk, since the <fi n form a complete set: 



(14) 



N g 

V" (j>n(xk)(fin{Xi) = —■ 

— ' Wi 
n—1 



It should be underlined that the m(x) defined by Eq. (|1£ 
are highly localized in the vicinity of the quadrature 
points but are not true delta- functions, since they are 
finite for x ^ Xk- 

Since the coordinate eigenfunctions Ui(x) are defined 
as continuous functions of x, it is possible to differen- 
tiate and tabulate them using the known properties of 
the (fin(x). The crucial step, of course, is to be able to 
find the set of basis functions and the related generalized 
quadrature which justify Eqs. (11,12]). 

While no conditions have thus far been placed on the 
set <fi n (x), other than orthonormality and completeness, 
in practice the requirement that the Ui(x) satisfy Eq. (|^) 
and Eqs. (0,0) limits the basis functions to the classical 
orthogonal polynomials. These are defined by a three 
term recursion relation of the form 

(ij(fij(x) = (x - aj)<f>j-i(x) - 0j-!<t>j-2(x), (15) 

with the properties 
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4>o(x) = const. 



(16) 



where the inner product is defined in Eq. (|T^). Since 
Eq. (15) may be interpreted as the definition of the coor- 
dinate operator x in the polynomial basis, the orthonor- 
mal eigenvectors of this tridiagonal matrix must diago- 
nalizc the coordinate operator pjfl . It is therefore not a 
coincidence that the associated eigenvalues are the gener- 
alized Gauss quadrature points associated with the orig- 
inal basis 4> n (x). For this set of N functions there is an 
associated generalized Gauss quadrature consisting of N 
points and weights which assure us that Eqs. (pj] , p!^ are 
satisfied exactly. 

There are as many ways to define the Ui(x) as there 
are classical orthogonal polynomials. The most natu- 
ral choice, however, is to use the Lagrange interpolat- 
ing functions, which explicitly diagonalize Eq. (|l5|). The 
Lagrange polynomials may be defined at the associated 
Gauss quadrature points by: 



JV 



Vi(x) = Yl 



fc=l 



Xk 



(17) 



where the prime denotes exclusion of the point Xi in the 
product. The Gauss-Legendre quadrature points Xk are 
defined on the interval [-1:1]; with the associated weights, 
one obtains the eigenfunctions: 



m(x) = Vi{x)/y/Wi- 



(18) 



These polynomials have the "delta-function property" by 
construction, and are normalized so that they form an 
orthonormal set under the generalized TV-point Gaussian 
quadrature: 



(ui | uj) = 6i,p 



(19) 



This definition ensures that there is an equality between 
exact integration and integration by quadrature for any 
integrand which may be represented as a polynomial of 
degree (2N — 1) or smaller. In the present work, both 
Lagrange functions as well as a DVR based on a Hermite 
polynomial basis are considered. 

Since the Ui(x) diagonalize the coordinate operator, 
the matrix element of any operator 0{x) which is a local 
function of x satisfies 



(«i I 0(x) | Uj ) = 5ijO(xi). 



(20) 



The DVR basis not only considerably simplifies the eval- 
uation of many matrix elements of the Hamiltonian, but 
leads to a sparse representation as well. For many large 
matrix problems the only practical methods of diagonal- 
ization or matrix inversion require the operation of the 
Hamiltonian matrix on some known vector. Sparsity is a 
key ingredient to performing this operation efficiently. 

While the result (^p| ) is an identity within a particular 
iVth-order Gauss quadrature, it is only exact when the 
product of Ui(x), uj(x) and the local operator 0(x) is a 



polynomial of degree 2 N — 1 or smaller Eq] . The present 
formalism, therefore, is particularly conducive to the so- 
lution of the GP equation in the limit of large No. In this 
case, the contribution of the kinetic term to the total en- 
ergy is small. The square of the condensate wavefunction 
may then be written in the so-called Thomas- Fermi (TF) 
approximation: 

V 2 (r) a [/ztf - \(x 2 + a 2 y 2 + /3 2 z 2 )] /4ir Vo , (21) 



where the normalization condition 
potential in the TF limit: 



i yields the chemical 



Mtf 



(15*) 2/5 = h (15ry) 2/5 . 



(22) 



in units of Twjq. In the TF limit, both the interaction and 
confining potentials appearing in the Schrodinger oper- 
ator L m) are evidently second-order polynomials. For 
a finite but large number of atoms, however, there is a 
small exponential tail near the boundary of the conden- 
sate cloud resulting from the finite kinetic energy Q [l8| . 
As shown in Sec. Ill, this exponential behavior may be ef- 



fectively captured even when using a comparatively low- 
order basis. 

For the present three-dimensional calculation, the con- 
densate and excited-state wavefunctions are expanded 
using Cartesian coordinates in a product basis of coor- 
dinate functions for each dimension; in the condensate 
case, for example, one writes: 



ip(r) = CjjkUj(x)uj(y)uk(z). 



(23) 



The components of the three-dimensional kinetic energy 
T x . VtZ = — 5 V 2 , y z has the following representation in the 
DVR product basis: 

(uiUjUk I T I uiu m u n ) 

= TijSj :m 8k,n + Tj >m Sij5k,n + Tk,nSi,lSj,m- (24) 

Thus, the Hamiltonian matrix separates into a sum of 
dense, one-dimensional kinetic energy matrices times 
Kronecker delta functions in the remaining variables, plus 
a purely diagonal matrix associated with the potential 
terms. While multidimensional expansions of this kind 
(^H) exist in other coordinate systems [^9|-^2) , the choice 
of a product basis in Cartesian coordinates yields a sepa- 
rable kinetic energy operator which ensures a sparse ma- 
trix representation of the Hamiltonian in the multidimen- 
sional DVR space. The expansion coefficients cy^ found 
by diagonalization are proportional to the values of the 
wavefunction at the appropriate quadrature points. 



B. Iteration and the DIIS Method 

The central numerical difficulty in the solution of the 
GP equation is the non-linearity associated with the two- 
body interactions. Since the magnitude of the Hartree 
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potential is proportional to the number of atoms in the 
condensate, the self-consistent solution when Nq is large 
can become problematic. Many calculations reported to 
date have "inverted" the search; that is, the value of /x is 
fixed while the wavefunction and associated value of Nq 
are obtained which satisfy the appropriate boundary and 
normalization conditions. Such a procedure is straight- 
forward in one dimension (though somewhat less so in 
two) |7j , since a direct numerical integration of the GP 
equation can be implemented. In three dimensions, a 
root search procedure may be devised to accomplish the 
same thing, but the calculation would be extremely time- 
consuming. With the method of steepest descents, the 
solution of the three-dimensional GP equation could be 
obtained by imaginary-time propagation. Though this 
approach has been shown to yield accurate results for 
large numbers of atoms Pq-|38|], it is not obvious that 
time-dependent techniques should be most suitable for 
the solution of an eigenproblem. 

An elegant and (as shown below) inexpensive method 
for the solution of the GP equation is to fix N and to 
self-consistently solve for fi and V'( r )- Unfortunately, a 
direct iteration of Eq. (|l|) does not usually converge to the 
global energy minimum, even for relatively small conden- 
sates (N ~ 10 3 , -q ~ 1). Initially, we attempted several 
simple schemes in order to improve convergence. One ap- 
proach was to use the TF approximation (Ell) to begin the 
iteration sequence. This would seem to be an excellent 
starting point when the condensate is large and the TF 
solution is a reasonable approximation. In practice, we 
have found that the poor behavior of the TF wavefunc- 
tion near the condensate boundary did not, in general, 
make this a viable procedure. Even correcting the TF 
result using a boundary-layer perturbation approach |4S|] 
did not seem to greatly improve convergence. 

Evidently, a robust numerical process is required which 
rapidly damps out the errors as the iteration sequence 
proceeds. One rather crude approach is to use a linear 
combination of the solutions from the (i — l) th and i th 
steps to initiate the iteration process for the (i + 1) 
step f55|]54|| ; the proper linear combination may be found 
by numerical experimentation. The approach used by 
the present authors, known as Direct Inversion in the 
Iterative Subspace (DIIS) p5[, is a more systematic ver- 
sion of the above that uses the information from all of 
the previous iterations. The DIIS method is well known 
in the quantum chemistry community, where it is used 
to accelerate the convergence of self-consistent field cal- 
culations. As discussed below, for small to intermediate 
numbers of atoms (10 3 < N < 10 6 ), DIIS is found to 
be at least as computationally efficient as the method of 
steepest descents. 

DIIS begins by defining an error which characterizes 
the convergence of the iteration process. While the error 
e may be defined in any number of ways, one that appears 
to be particularly stable numerically is the commutator: 



where p(r,r') is the density matrix Evidently, e 

vanishes at convergence. Note that both the density ma- 
trix and the error are matrices having S 2 elements, where 
S is the dimension of the Hilbcrt space. The procedure is 
to expand the current expression for the Hamiltonian and 
error as a linear combination of the m previous values 



i=l 
m 

i=i 

subject to the constraint that 

m 

$3«n = i. 



(26) 
(27) 

(28) 



i=i 



Minimizing the squared norm of i m+1 with the con- 
straint (Eq), one obtains m + 1 linear equations: 



X! '''' • A °- 



where 



B 



(29) 



(30) 



and j — 1,2, ... ,m. The Lagrange multiplier A enforcing 
the constraint yields the squared norm of the error. 

The practical implementation of the DIIS algorithm, 
especially for large Hamiltonian matrices, deserves some 
additional comment. Although it appears to be neces- 
sary to store the Hamiltonian matrix for each iteration, 
this is not the case. The only part of the Hamiltonian 
matrix that varies from iteration to iteration is the non- 
linear potential, and this matrix is diagonal in the DVR. 
Inserting Eq. (|30| ) into (j2^) , the matrix elements Sj j are 
simply given by: 



B itj = 2(cpi | ipj)(ipi | ipj) - 2(cpi | ipj){ipi | tpj) 



(31) 



(25) 



where tpi(r) = H l ipi(r); it is important to note that H l is 
the Hamiltonian matrix constructed from -0i_i(r), whose 
minimal eigenvector is ipi{ r )- O n each iteration, there- 
fore, it is only necessary to store the solution vector as 
well as the vector which represents the operation of the 
Hamiltonian matrix on the solution vector. Since these 
objects are both of dimension S rather than S 2 , they may 
easily be stored in central memory even on current work- 
stations as long as the number of DIIS iterations does 
not become too large. 

To summarize, only the vectors ^(r) and ip(r) need to 
be stored at each iteration step. From these, the scalar 
products ( |3l| ) required to set up the set of linear equa- 
tions ( p8| ) and ( p9| ) may be generated. Once the m un- 
knowns a,i as well as A have been found, a new nonlinear 
potential may be constructed from the past guesses, 
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Vg* 1 = Y,«iV&, (32) 

i=i 

and the next self-consistent field cycle initiated. Further- 
more, by storing and not destroying all previous elements 
of Bi j , only an additional row and column need be com- 
puted at each step of the iteration. Since the dimension 
of B is assumed to be small, the additional storage should 
not be a problem on most workstations. 

With the TF expression for the condensate wavefunc- 
tion as the initial guess, the DIIS algorithm yields a con- 
vergent solution for significantly larger numbers of atoms 
(No ~ 10 5 ) than was possible with a direct iterative pro- 
cedure. As No continues to grow, however, the number 
of DIIS iterations to achieve convergence increases until 
eventually the process fails. It appears that small errors 
in the intermediate solution of the GP equation, partic- 
ularly in the 'tail' region close to the condensate surface, 
are amplified by the nonlinearity of the potential during 
the iteration process. Combining DIIS with a more real- 
istic initial guess, we have been able to achieve conver- 
gence for the large number of atoms (No ~ 10 6 ) relevant 
experimentally. Two schemes used to improve the start- 
ing approximation for the Hartree potential have proved 
particularly valuable. By slowly increasing the number 
of atoms, the converged solution for a slightly smaller 
condensate is an excellent choice. Alternatively, one may 
solve a modified GP equation with an exaggerated kinetic 
energy contribution, then slowly ramp down the degree 
of exaggeration (this technique is similar in spirit to sim- 
ulated annealing). A simple modification of the original 
DIIS procedure was also found to be quite helpful. The 
DIIS procedure is broken into cycles with some maxi- 
mum number of iterations allowed per cycle. At the end 
of each cycle, the DIIS procedure is restarted using the 
best available solution. By numerical experimentation it 
was found that once the root mean square error (or al- 
ternatively \/A) is reduced to 10~ 2 — 10~ 3 , it becomes 
possible to restart DIIS and to rapidly reduce errors to 
10~ 6 — 10~ 7 with very few additional iterations. 

C. Interpolation of the Wavefunction 

As discussed above, in the limit of very large No where 
the TF theory is believed to be valid, both the interac- 
tion and confining potentials are well-approximated by 
low-order polynomials in x, y, z. Since an AT-point Gauss 
quadrature is able to integrate a (2N — l)th-order func- 
tion exactly, it should be possible to capture the essential 
behavior of the condensate wavefunction (squared) using 
only a very small number of quadrature points. One 
is therefore left with a somewhat surprising conclusion: 
when the number of atoms becomes very large, the re- 
sulting self-consistency problem should simplify consid- 
erably. A very coarse DVR grid has two obvious advan- 
tages. The reduced dimension of the Hamiltonian matrix 



would accelerate the solution of the eigenproblem at each 
DIIS iteration. In addition, the fewer degrees of freedom 
for spatial variations should enable the self-consistent so- 
lution for the condensate to converge with fewer DIIS 
iterations. The crucial unknown is whether the exponen- 
tial tail at the condensate boundary, due to the finite 
kinetic energy contribution, varies sufficiently rapidly to 
invalidate the low-order approximation applicable in the 
TF limit. 

Suppose that the gross features of the condensate 
wavefunction could be found using a low-order Gauss 
quadrature, corresponding to only a few DVR points in 
each spatial direction. The self-consistent solution ob- 
tained using such a 'coarse grid' would then make an 
excellent initial guess for a more accurate 'fine grid' cal- 
culation, if the interpolation between grids could be im- 
plemented successfully. By increasing the number of 
DVR points in the mesh once or perhaps a few times, it 
should be possible to rapidly converge the solution of the 
GP equation for virtually arbitrary numbers of atoms. 
The interpolated wavefunction on the quadrature points 
(xi,y m , z n ) associated with the fine grid is approximately 

ip(xi,y m ,z n ) w CijkUi(xi)uj(y m )uk(z n ), (33) 

ijk 

where the sum is taken over coarse-grid expansion coef- 
ficients and coordinate eigcnfunctions. Since the DVR 
points corresponding to the two different Gauss quadra- 
tures are not generally coincident, the values of the 
coarse-grid coordinate basis functions at the fine grid 
points must be obtained explicitly using either Eq. <\12\ ) 
or @. 

D. Extraction of Eigenvalues and Eigenvectors 

The direct extraction of all the eigenvalues and eigen- 
vectors of Eqs. (Q-§) becomes computationally imprac- 
tical for very large matrices. In three-dimensional sys- 
tems, a sizeable number of basis functions (and therefore 
quadrature points) are usually required in order to ade- 
quately represent the wavefunctions over all space. This 
is a particularly important consideration for the excita- 
tions, whose rapid spatial variations can only be captured 
by high-order polynomials. As a result, the dimension of 
the matrices can become very large, even after block- 
diagonalizing according to parity. A full diagonalization 
of the eigenproblem, using standard and widely-available 
routines based on variants of the Givens-Householdcr 
method, places considerable demands on both storage 
and cpu-time. A full diagonalization is not necessary 
since only the lowest eigenvalue (fj,) of the GP equation 
is required, and it is sufficient to determine merely the 
lowest-lying excitations of the Bogoliubov equations at 
zero temperature. 

When the matrices are too large to fit in memory, it- 
erative methods must be used in order to extract the rel- 
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evant low-lying states. These techniques require the fre- 
quent operation of the Hamiltonian on a vector, and are 
practical only if there are relatively few non-zero matrix 
elements. Indeed, a sparse representation of the Hamilto- 
nian is a crucial feature of the DVR and the separability 
of the kinetic energy operator. In addition, standard iter- 
ative procedures return reasonable approximations only 
to the extremal eigenvalues, in the spirit of a Rayleigh- 
Ritz variational method. The GP equation ([[]) yields 
a positive-definite, real, symmetric, and sparse matrix 
whose lowest eigenvalue may be found using variants of 
the well-known Lanczos |i5| or Davidson |^6| algorithms. 

The Davidson method, which is widely used by quan- 
tum chemists, develops a small orthonormal subspace of 
vectors which is adequate to describe the eigenpairs of 
physical interest. The Hamiltonian is projected into this 
p-dimensional subspace and the Rayleigh-Ritz variational 
principle ensures that an upper bound to the lowest n 
eigenvalues is obtained from the process. Since the com- 
putational procedure in the application of the Davidson 
algorithm may be unfamiliar to some readers, it is sum- 
marized below: 

1. Choose an initial set of m orthonormal trial vectors 
bj , where m > n. 

2. Calculate hj = Hhi, the effect of the Hamiltonian 
matrix on these vectors. 

3. Calculate the Hamiltonian matrix (b{|hj) from the 
information generated in steps 1 and 2. 

4. Solve the resultant small eigenvalue problem for 
the current values of the wavefunctions, \I>£ UI and 
eigenvalues, Eg UI by a direct method. 

5. Calculate the residuals, 

r q = (E™ - ff) * c g UI ' (34) 

and test for convergence. If all roots are converged 
stop. If not select the unconverged residuals for 
further improvement. 

6. Using the set of residuals solve the equation, 

(fi - E™) * q = r q (35) 

where H is some approximation to the exact H. 

7. Schmidt orthonormalize the <J> 9 to the and then 
append them to the b^ to enlarge the vector space. 

8. Calculate = Hhi, for the appended vectors and 
the additional matrix elements needed to border 
(bilk,-). 

9. Return to step 4. 



In the original Davidson algorithm, which was de- 
signed for diagonally dominant matrices, the approxi- 
mate Hamiltonian used in step 6 was the diagonal. Then 
the calculation of the solution to Eq. (|35| ) is trivial. For 
diagonally dominant matrices, this Jacobi precondition- 
ing method produces a quickly converging set of new vec- 
tors, and the entire calculation is dominated by step 2 in 
the sequence above. Although the DVR has the desir- 
able feature that it produces a sparse representation, the 
matrix is unfortunately far from being diagonally domi- 
nant. Without an alternative preconditioning technique, 
the calculation can require several (~ S) full cycles which 
becomes computationally prohibitive. 

Our approach to the preconditioning problem is based 
on using a separable approximation to H, i.e. a Hamilto- 
nian which can be written as a sum of commuting oper- 
ators in the three-dimensional space. For simplicity, we 
have chosen to use the bare trap Hamiltonian, but other 
more complicated choices can be made, subject to the 
separability requirement. Eq. ( |35| ) then becomes: 

(nl + nl + nl-E™)<f> q = r q . (36) 

In order to solve this equation, we transform to the 
basis set which diagonalizes the separable Hamiltonian, 
recognizing that the transformation can be written as 
a product of three one-dimensional unitary transforma- 
tions. Thus, the residual vector on the right hand side 
of Eq. ( |34j ) is transformed from the DVR to the diagonal 
representation, an operation which scales as the size of 
the three-dimensional basis set times the one-dimensional 
basis set. This linear scaling with basis set size is a 
consequence of the separable form of the approximate 
Hamiltonian, yet preserves the most desirable features 
of the DVR. Once the residual has been transformed to 
the diagonal representation, Eq. ( |35"1) is trivially solved 
and then the solution vector is transformed back to the 
DVR representation. This procedure is analogous to us- 
ing fast Fourier transforms (FFT) to solve multidimen- 
sional PDE's. The FFT is more efficient (but less gen- 
eral) due to special character of the Fourier transforma- 
tion but the separability is a key factor to its proper- 
ties. By using this separable preconditioning, we have 
been able to reduce the number of Davidson iterations 
to manageable size. In addition, the vectors from a pre- 
vious calculation for a smaller condensate are found to 
be excellent trial vectors to initiate the calculation. The 
result is a robust procedure which is quite efficient for 
large matrices. 

The matrix associated with the Bogoliubov equa- 
tions (d||), in contrast, is non-symmetric and has real 
eigenvalues occurring in positive-negative pairs. This 
matrix must be addressed using routines such as those 
based on an Arnoldi or a non-symmetric Davidson ap- 
proach. Since the low-lying excitations of the Bogoliubov 
equations are those closest in magnitude to zero (rela- 
tive to the chemical potential) and are therefore the least 
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extremal, it is most convenient to decouple the original 
equations using the following linear combinations: 

/„(r) = u n (r) + v n {r); (37) 
g n {r) = u n (r) - v n (r). (38) 

The resulting eigenproblem then becomes 

(L - ii + 2Vh) (t - fx) /„(r) = e 2 n f n (r), (39) 

where eo — corresponds to the ground-state energy 
relative to the chemical potential. Once the functions 
/„ are found, the left eigenvectors g n may be obtained 
directly by utilizing 

g n (r) = e- 1 (l - fi) /„(r) (40) 

subject to the normalization condition (f\g) — 1, ob- 
tained from Eq. (||) . While the resulting eigenproblem re- 
mains non-symmetric (since the condensate density does 
not commute with the operator L), the dimension of the 
matrix has been reduced by a factor of two. Further- 
more, the eigenvalues are positive-semidefinite and 
real ||. In practice, the composite operator on the left 
side of Eq. ( p9|) is never determined explicitly; doing so at 
each DIIS iteration not only would be time-consuming, 
but also would lead to a dense matrix. Rather, at each 
Arnoldi/Davidson iteration the vector is multiplied in 
turn by the two distinct sparse operators. 

III. RESULTS AND DISCUSSION 

The techniques discussed in the previous section are 
particularly useful under two circumstances: when the 
potential is completely anisotropic, and when the number 
of condensate atoms is very large. For illustrative pur- 
poses, results arc presented for Bose-condensed sodium 
atoms confined in traps with three different geometries, 
all of which have been realized experimentally. The 
completely anisotropic trap considered here is a TOP 
trap with angular frequencies in the natural ratio |57]] 
(uj x ,uj y ,uj z ) = u)q(1, \/2, 2) where = 3547T rad/s; BEC 
in such a system has been recently observed by Phillips et 
al. |i8| . Other geometries considered are the cigar-shaped 
Ioffe-Pritchard trap of W. Ketterle |59| with cylindri- 
cal symmetry and (u x , uj y , ui z ) « ^(l, 13.585, 13.585), 
where ujq w 33.867T rad/s, and the approximately spher- 
ical '4D' potential of L. Hau [Q with w 87 rad/s. All 
the calculations assume an s-wave scattering length for 
sodium of a = (52 ± 5)ao, where ao is the Bohr radius 




A. Condensate 

The solution of the GP equation ([!]) on the coarse 
and fine grids employed 60 and 180 harmonic oscillator 



basis functions (Hermite polynomials with a Gaussian 
prefactor) in each spatial direction, respectively. Coor- 
dinates were rescaled through {x,y, z} — > {x, y/a, z//3}; 
the non-linear coefficient becomes g — Airr] — Ana(3riQ 
in order to ensure the proper normalization of the con- 
densate density (0). With this choice, the condensate 
cloud becomes almost spherical at large particle numbers, 
with an approximate radius given by the TF expression 
R TF = (I577) 1 / 5 , in units of the characteristic trap length 
da = ^/fi/Mujo- DVR points significantly beyond the TF 
radius were ignored. The ground state was assumed to 
have totally even parity, so the GP equation was solved 
in a single octant. The basis functions were cigenfunc- 
tions of a trap with frequencies reduced from their actual 
values by a factor of 20, in order to decrease the density 
of points. The resulting Hilbert spaces for the coarse 
and fine grids then have dimension 3 566 and 18 685 re- 
spectively, for all geometries and condensate numbers; 
both sizes are considerably reduced from the impractical 
values of 216 000 and 5 832 000 one naively might have 
obtained. 

In Table |, the chemical potentials obtained numer- 
ically using the coarse grid are given as a function of 
the number of atoms A*o = 2 q in the condensate for 
the spherical, cylindrical, and anisotropic geometries de- 
scribed above. The TF values, from Eq. (|22|), are in- 
cluded for comparison. The number of DIIS iterations 
required to yield a convergent solution to the GP equa- 
tion increases with the number of atoms; in the cylindri- 
cal case, DIIS failed to converge for N = 2 19 = 524 288 
and greater. Choosing a finer grid always reduces the 
number of DIIS iterations; although increasing the num- 
ber of points yields additional degrees of freedom for 
the Hartree potential, the variations of the condensate 
density (particularly in the surface region) are better 
captured. In the cylindrical case, which has a cigar 
shape, the coarse grid has too few DVR points in the 
two strongly-confining directions to adequately capture 
the behavior of the condensate at large N . For q > 18, 
therefore, initial grids were chosen to contain a larger 
number of points not exceeding 10 4 . Since more points 
imply a larger eigenproblem at each iteration, practical 
calculations require as coarse a grid as possible. 

The convergence of the chemical potentials with the 
number of condensate atoms is shown in Fig. [|. Both the 
cylindrical and spherical traps give rise to relatively slow 
convergence of the chemical potential to the TF value 
compared with the anisotropic case. This trend is ex- 
pected for the spherical trap, where the confinement is 
extremely weak. In the cylindrical case the condensate is 
strongly confined, but only in the radial direction where 
the motion of the atoms is more or less frozen. Since 
the axial kinetic energy can be large, however, the cylin- 
drical trap is effectively loose and one-dimensional. In 
contrast, the fully anisotropic trap considered here is rel- 
atively tight in all directions, and the chemical potential 
converges to the TF value more rapidly. 
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FIG. 1. The percent difference between the numerical and 

TF values for the chemical potential A/j, = 1 — /^tf / Mcxact are 

shown as a function of the number of atoms and the universal 

scaling parameter rj = a/3(Noa/do) (which is proportional to 

the non-linear coupling constant) for all three trap geometries. 

The data for the spherical (circles), cylindrical (diamonds), 

and anisotropic (squares) cases are identical to those given 

in Table M. The solid, dashed, and dotted lines are fits to 

{S,C,A} /( ,{S,C,A} 



the data points using the expression Aji ~ 7 
for the spherical (7 s = 1.5), cylindrical (7 



c 



170), and 

anisotropic (7^ = 4.0) cases, respectively. It is important to 
note that in the inset, the values of A^t for the cylindrical data 
are reduced by a factor of 10 in order to facilitate comparison. 

In the TF limit, most relevant quantities, such as the 
mean condensate radius i?TF = (15t7) 1/,5 c?o and the chem- 
ical potential /iTF = are functions of the 
universal scaling parameter rj = a(3(Noa/do). Simi- 
larly, the first-order correction to the TF chemical po- 
tentialuxF, taking into account the average kinetic en- 
ergy [|46|,[l7[ and contributions from the potential ener- 



gies [1481, is proportional to -R^r ~ Mtf- The ^ s t° the 
numerical data of A/i = (1 — /itf/m) ~ 7/mtf, shown 
in Fig. [I], are in reasonable agreement with this behavior 
(except for the cylindrical case for very small numbers 
of atoms iVo < 10 4 ). One might naively expect, there- 
fore, that the convergence of the chemical potential to 
the TF limit as a function of rj (which is proportional 
to the coefficient of the non-linear term in the GP equa- 
tion) would be independent of trap geometry. As shown 
in Fig. [l], this is in fact not the case. The data for the 
cylindrical case converge far more slowly than those for 
the other geometries, though the magnitudes of the TF 
chemical potentials for a given value of rj are identical 
(note that the values of A/z shown in the figure for the 
cylindrical case are reduced from their actual values by 
an order of magnitude in order to facilitate comparison). 
The kinetic energy contribution to the chemical potential 
(and therefore to the total energy) is evidently strongly 
dependent on trap geometry. Thus, a large non- linear 
term in the GP equation does not necessarily mean that 
the TF approximation adequately represents the system. 
In order to verify the accuracy of the solution obtained 



on the coarse grid, the condensate wavefunction was in- 
terpolated onto a 18 685-point grid derived from a DVR 
basi s wit h 180 polynomials in each direction (refer to 
Sec. II C for details on the interpolation technique). The 
solution of the GP equation was again converged on the 
fine grid. For small to intermediate values of the nonlin- 
ear coupling (77 < 3000), a single interpolation between 
the two meshes was sufficient to yield a solution on the 
fine grid with only a few (~ 10) additional iterations. As 
77 increased, two or more successive interpolations and 
re-convergences were generally required before the solu- 
tion on the finest grid could be obtained. In all cases, the 
values of the chemical potential at the two extremes were 
found to be identical to at least three decimal places. 

The condensate wavefunctions obtained on the coarse 
and fine grids are compared in Fig. |2| for the case of 
Nq = 2 19 = 524 288 atoms in the completely anisotropic 
trap. The coarse-grid condensate profile along a given 
axis appears rather crude, particularly in the surface re- 
gion where the wavefunction varies rapidly; for a given 
basis, convergence is only obtained at the DVR points. 
Nevertheless, this wavefunction interpolated is virtually 
indistinguishable from the converged solution on the fine 
mesh. As rj increases, however, the length of the tail at 
the surface of the condensate shortens. Eventually, the 
coarse grid will have too few points in this crucial region 
to adequately capture the rapid variations. In this case, 
a large jump in mesh size results in an interpolated wave- 
function that more poorly represents the self-consistent 
result. 
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1.500 


14.085 


2.207 


10 


1.825 


(1.119) 


17.384 (9.393) 


3.572 


(2.824) 


11 


2.065 


(1.477) 


19.392 (12.395) 


4.345 


(3.726) 


12 


2.435 


(1.949) 


22.359 (16.355) 


5.425 


(4.917) 


13 


2.970 


(2.571) 


26.620 (21.580) 


6.904 


(6.488) 


14 


3.719 


(3.393) 


32.682 (28.475) 


8.900 


(8.560) 


15 


4.743 


(4.477) 


41.055 (37.573) 


11.572 


(11.296) 


16 


6.124 


(5.907) 


52.433 (49.578) 


15.128 


(14.904) 


17 


7.970 


(7.795) 


67.750 (65.418) 


19.847 


(19.667) 


18 


10.427 


(10.285) 


88.228 (86.320) 


26.096 


(25.950) 


19 


13.685 


(13.571) 


115.453 (113.900) 


34.358 


(34.241) 


20 


17.999 


(17.907) 


151.545 (150.29) 


45.275 


(45.182) 


21 


23.702 


(23.628) 


199.313 (198.31) 


59.693 


(59.618) 


TABLE I. 


The chemical potentials jjS s ' C ' A ^ corresponding 



to spherical, cylindrical, and anisotropic geometries respec- 
tively, are given in units of Tilo\ S,C ' A ^ for various numbers of 
condensate atoms No = 2 q . The TF values, from Eq. (p2[), 
are given in parentheses. All results are converged to three 
decimal places and were obtained using the coarse grid with 
at least 60 3 basis functions and 3 566 DVR points. 
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{x,y,z} (units of d ) 
FIG. 2. The condensate wavefunction for 
No = 2 19 = 524 288 atoms in the fully anisotropic trap, nor- 
malized to unity, is shown as simultaneous projections along 
the positive x (rightmost curves), y (middle curves), and z 
(leftmost curves) axes. The dashed lines correspond to the TF 
approximation. The numerical results obtained on the coarse 
and fine grids are shown as dotted (offset 0.002) and solid 
(offset 0.004) lines, respectively. The interpolated coarse-grid 
and converged fine-grid wavefunctions exactly coincide. 



B. Excitations 

The excitations of a condensate in a fully anisotropic 
harmonic trap have been completely classified [pl|, and 
have been explicitly obtained in the low-density |6Sj] and 
TF |33],[34| limits. The states are polynomials of order 
N = I + m + n and are labeled by the total parity 
P = (— iy+ m + n ^ where the quantum numbers (l,m,n) 
represent the order of the polynomials along the (x, y, z)- 
directions in the non-interacting limit. In the strongly- 
interacting (or hydrodynamic) regime, there are four 
odd and four even parity low-lying modes with energies 
e = y/l + a 2 m + P 2 n in units of hu>o, where (I, in, n) can 
be either or 1. The ground state, relative to the chem- 
ical potential, has quantum numbers (0,0,0). The only 
states with N = 1 are the odd-parity dipole modes, where 
the center of mass oscillates with the three trap frequen- 
cies. For N = 2, there are six quadrupole oscillations 
with even parity and stationary center of mass. Three of 
these have energies given by the expression above, with 
(l,m,n) — (1,1,0) and cyclic permutations. The other 
three are the solutions of the secular equation: 



3-e 2 
1 
1 



1 



2 /a 2 



1 



V/9 2 



= 0. 



(41) 



For (a,f3) = (V2, 2), the geometry considered here, one 
readily obtains e = a/8± 4y/2 and y/E. In the fully 
anisotropic case, all the excitation energies (except those 
for the odd-parity dipole modes) decrease with Nq 163] . 
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q (JV„=2') 

FIG. 3. The low-lying excitations of a condensate in the 
completely anisotropic trap are given in trap units Tiujq as a 
function of q, where the number of atoms No = 2 q . The cir- 
cles correspond to numerical results; horizontal dashed lines 
are the predictions of the TF theory. The data points along 
zero are the ground-state energies relative to the chemical po- 
tentials listed in Table |l| The long-dashed modes are labeled 
by (Imn), where I, m, and n are the quantum numbers of 
the non-interacting harmonic oscillator states along x, y, and 
2, respectively. The unlabeled short-dashed excitations with 
energies just below and above 2hu>o are higher-order modes 
with odd parity along x and y, respectively; the two at higher 
energies have totally even parity (lower) and odd parity in 
both x and y (upper). 

The low-lying excitation energies for a completely 
anisotropic trap have been computed numerically, and 
are shown in Fig. ||. All the calculations were obtained 
using the 3 566-point grid. While small finite-size and 
coarse-graining effects are present, as evidenced by the 
small fluctuations in the ground-state energy e = 0, the 
results closely match the TF predictions described above. 
As expected, the odd-parity dipole modes are indepen- 
dent of the number of atoms in the trap. All the other 
frequencies depend strongly on No, decreasing from their 
non-interacting values in agreement with perturbative 
calculations |62j. 

While the frequencies of all the lowest-lying modes at- 
tain their large-number values by Nq ~ 10 6 in the rela- 
tively strong anisotropic trap, the convergence to the hy- 
drodynamic limit slows as the quantum numbers increase 
and the confinement is weakened. In Fig. ||, the low-lying 
excitations e n i of the loose spherical trap are shown as a 
function of the number of atoms in the condensate, up to 
N = 2 22 w 4.2 x 10 6 . In the hydrodynamic limit, the en- 
ergies are given by e n i = y/l + 2n(n + I + 3/2) ||. The 
lowest number-dependent mode £02 agrees with its hydro- 
dynamic value V2 to less than a percent by Ao « 10 6 . 
The numerical value of e Q4 , in contrast, differs from its 
limiting value of 2 by approximately 2% even when the 
number of atoms is as large as 4 x 10 6 atoms. Evidently, 
the magnitude of the excitation energy, relative to the 
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ground state //, does not alone provide a sufficient in- 
dication of its convergence to the hydrodynamic limit. 
A similar number-dependence for the higher-lying exci- 
tations of a cylindrically-symmctric condensate has also 
been recently obtained p3j. 
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FIG. 4. Selected low-lying excitations e n i of the spheri- 
cal condensate are given in trap units Tlujo as a function of 
q, where the number of atoms No = 2 9 . The filled and 
open circles correspond to numerical results for n = and 
n = 1, respectively, while squares represent £20; horizon- 
tal dashed (n — 0), long-dashed (n = 1), and dot-dashed 
(n = 2) lines are the results obtained in the hydrodynamic 
limit e nl = y/l + 2n(n + I + 3/2). 



to employ more sophisticated techniques, including num- 
ber or kinetic energy ramping and multigrid interpola- 
tion. In general, DIIS becomes more expensive as the 
number of atoms increases; it is conceivable that an al- 
ternative method such as imaginary-time propagation be- 
comes more efficient than DIIS in the regime Nq 10 6 . 

The convergence of the chemical potential and the low- 
lying collective excitations is investigated as a function 
of trap geometry. The chemical potential is found to 
approach its Thomas-Fermi value more slowly as the 
confinement weakens or the degree of anisotropy be- 
comes more pronounced; the convergence does not scale 
universally with the magnitude of the non-linear coeffi- 
cient. The excitations of a completely anisotropic con- 
densate have been calculated numerically, and for large 
numbers of atoms agree with the Thomas-Fermi predic- 
tions (6^^J] . For a very weak spherical trap, the collec- 
tive frequencies converge to their hydrodynamic values 
more slowly. 
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IV. CONCLUSIONS 



A numerical procedure is introduced for the investi- 
gation of an interacting Bose gas at zero temperature 
confined in a completely anisotropic trap. The central 
feature of the technique is the use of the discrete variable 
representation (DVR) as the primary basis for the cal- 
culations. The DVR combines the best features of grid 
and basis-set techniques. All local operators are diago- 
nal, so the evaluation of interaction matrix elements be- 
comes trivial. While the kinetic energy has a more dense 
representation, it may be evaluated analytically to high 
accuracy by exploiting the underlying polynomial basis 
used to define the DVR. Furthermore, the kinetic energy 
needs to be evaluated only once. 

In the present method the condensate density is deter- 
mined self-consistently; for fully three-dimensional sys- 
tems, this approach is considerably more efficient than 
conventional root-search algorithms. At each iteration, 
the ground-state wavefunction is obtained by iteratively 
diagonalizing the sparse GP Hamiltonian, using either a 
Lanczos or Davidson method. Convergence of the self- 
consistent solution to the GP equation is substantially 
hastened by employing direct inversion in the iterated 
subspace (DIIS). As the non-linear coefficient of the GP 
equation becomes very large, it often becomes necessary 
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